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ABSTRACT 

The applicability of the highly idealized secondary infall model to ‘realistic’ initial 
conditions is investigated. The collapse of proto-halos seeded by 3cr density perturbations 
to an Einstein-de Sitter universe is studied here for a variety of scale-free power spectra 
with spectral indices ranging from n = 1 to —2. Initial conditions are set by the con¬ 
strained realization algorithm and the dynamical evolution is calculated both analytically 
and numerically. The analytic calculation is based on the simple secondary infall model 
where spherical symmetry is assumed. A full numerical simulation is performed by a Tree 
N-body code where no symmetry is assumed. A hybrid calculation has been performed by 
using a monopole term code, where no symmetry is imposed on the particles but the force 
is approximated by the monopole term only. The main purpose of using such code is to 
suppress off-center mergers. In all cases studied here the rotation curves calculated by the 
two numerical codes are in agreement over most of the mass of the halos, excluding the 
very inner region, and these are compared with the analytically calculated ones. 

The main result obtained here, reinforces the foundings of many N-body experements, 
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is that the collapse proceeds ’gently’ and not via violent relaxation. There is a strong 
correlation of the hnal energy of individnal particles with the initial one. In particnlar 
we hnd a preservation of the ranking of particles according to their binding energy. In 
cases where the analytic model predicts non-increasing rotation cnrves its predictions are 
conhrmed by the simnlations. Otherwise, sensitive dependence on initial conditions is 
fonnd and the analytic model fails completely. In the cosmological context power spectra 
with n > —1 yields (in the mean) non-increasing rotation cnrves, and in snch cases the 
secondary infall model is expected to be a nsefnl tool in calculating the hnal virialized 
structure of collapsing halos. 


I. INTRODUCTION 

In a cosmological model in which structure forms via gravitational instability, the 
collapse of proto-structures onto local density maxima of the primordial perturbation held 
plays a major role in structure formation. Although the hnal (non-linear) density maxima 
do not necessarily emerge out from initial (linear) density maxima (c/., Hohman 1986, and 
1989, Katz et a/. 1993, Zaroubi and Hohman 1993b and Bertschinger & Jain 1994) there 
is a good evidence to believe a strong connection between initial high maxima and hnal 
density maxima (van der Weygaert & Babul 1994). 

The problem of collapse seeded by such maxima has been investigated from two dif¬ 
ferent points of view. One is the statistical distribution of the objects formed this way, 
namely the distribution of galaxies and clusters. This is related to the so-called biasing 
problem which has been extensively studied (Kaiser, 1984, Davis et al. 1985 Bardeen et 
al. 1986 hereafter BBKS). The other aspect of the collapse seeded by local peaks is the 
internal structure of the objects thus formed and its dependence on the statistical proper¬ 
ties of the primordial perturbation held. This problem had been hrst addressed by Gott & 
Gunn (1972) and was followed by the seminal paper of Gunn (1977, hereafter G77). The 
ideas and formalism suggested by G77 have led to considerable ehorts to understand the 
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dynamical and statistical aspects of the problem, nsing both analytical calcnlations and 
n um erical N-body simnlations to nnderstand this highly non-linear collapse process. It is 
the aim of the present paper to present N-body simnlations which have been designed to 
check and make detailed comparison with the predictions of the analytical calcnlations. 

Gnnn & Gott (1972) stndied the very idealized problem of the collapse onto a single 
feirac-like density pertnrbation to an otherwise nnpertnrbed flat Friedmann nniverse. This 
process has been described later on as the secondary infall, namely the collapse of less 
bonnd shells onto a pertnrbation that has already collapsed and virialized. By following 
the trajectories of spherical shells, and ignoring the process of shell crossing, the hnal 
virialized strnctnre was calcnlated. The major breakthrongh was done by G77 who realized 
that the spherical collapse onto the peak admits an adiabatic invariant. The case of a self¬ 
similar pertnrbation embedded in a flat Einstein-de Sitter nniverse allows an analytical 
calcnlation of the hnal virialized strnctnre which takes into acconnt shell crossing. The 
pioneering work of G77 was extended by Filmore and Goldreich (1984, hereafter FG) and 
Bertschinger (1985). These anthors presented a new nnmerical approach to the problem 
of self-similar collapse in a hat nniverse. These self-similar solntions provide an exact 
description of the collapse all the way to the hnal virial eqnilibrinm. FG farther extended 
the analytical approach of G77 and presented an asymptotic analytical solntion of the 
density strnctnre, the natnre of which is described below. A similar analytical calcnlation 
was presented also by Bertschinger and Watts (1988). 

The problem of collapse and in particnlar the formation of dark halos, hereafter this 
term describes bonnd objects made of collisionless matter, was stndied in great detail by 
Qninn, Salmon and Znrek (1986; QSZ), Znrek, Qninn and Salmon (1988), Frenk et al. 
(1988) and Grone, Evrard and Richstone (1994). These anthors addressed the problem 
by means of N-body simnlations, in which the fnll non-linear evolntion from typical initial 
conditions corresponding to diherent models was followed. The main emphasis of these 
papers is on the dynamical resolntion of forming objects and the stndy of their strnctnre. 
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Thus, unlike the analytical approach where very idealized physical systems have been 
considered, the N-body simulations enable the study of the evolution of typical objects 
whose structure is not self-similar and does not obey a high symmetry. These studies 
focused on hnding systematic trends in the forming structures, and in particular possible 
dependence on the initial conditions. 

In the standard model of cosmology structure emerges out of a primordial perturba¬ 
tion held. This held is assumed to constitute a random Gaussian held, whose statistical 
properties are dehned by its power spectrum. A hrst step to bridge the gap between 
the highly idealized analytical calculations and the numerical simulations was taken by 
Hohman and Shaham (1985, hereafter HS), who calculated the ensemble mean density 
prohle around high local density maxima. Assuming that structure around peaks of the 
primordial held is indeed given by the ensemble average, and that the dynamics is well 
approximated by the G77 dynamical model, HS calculated the virialized hnal structure 
expected in a variety of cosmological models. The analytical calculations of HS present a 
highly simplihed picture of the extremely asymmetrical and very non-linear actual physi¬ 
cal process. In spite of this high level of simplihcation, the predictions made by HS have 
been basically conhrmed by N-body simulations (QSZ, Frenk et al. 1988). Hohman (1988) 
made a detailed comparison of analytical predictions of the so-called secondary infall model 
(hereafter SIM) with the simulations of QSZ and Quinn and Zurek (1988) and found a 
good agreement in terms of the rotation curves of dark halos that are formed in the sim¬ 
ulations. Thus, the simple model successfully reproduces the virialized structure of halos. 
However, a closer inspection of the formation process of these halos reveals a picture that 
stands in complete disagreement with the SIM. Halos do not form by the infall of spherical 
shells, but rather by the coalescence and mergers of substructures (QSZ). Nothing of the 
symmetric and ordered process envisaged by HS is actually co nfi rmed by the numerical 
simulations, in which the collapse process of halos seems to be an extremely unordered 
and random process. In spite of these two very different formation pictures, the analytic 
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and numerical calculations yield a similar final structure. It is this seemingly paradoxical 
behavior that the present paper addresses, by designing and running N-body calculations 
so as to probe the basic ingredients of the collapse process of halos. The main aim of the 
paper is to hud the key to the success of the SIM, and at the same time to determine the 
limitations of this simple model. A major goal here is to dehne under what conditions 
the model can be applied to realistic initial conditions and to serve as a practical model 
for bound structure formation in an expanding universe. An entirely different approach to 
the problem has been adopted by Quinn and Zurek (1988) and Warren et al. (1991), who 
suggested that the basic mechanism is related to the mergers of subclumps rather than 
centralized collapse. The role of mergers will be studied here by the N-body simulations. 
The structure of the paper is as follows. In §II the main features of the analytical model 
are presented and the basic question which are to be solved by the numerical simulations 
are posed. This is followed by the description of the numerical simulations, in §III, and 
the initial conditions, in §IV. The results of the N-body simulations are analyzed in §V. 
The paper concludes with a discussion given in §VI. 

II. Secondary Infall 

The main essence of the SIM is that for a single spherical scale-free perturbation to an 
otherwise homogeneous flat universe, the hnal time averaged radius of a given Lagrangian 
shell scales with its initial radius. Thus, for a primordial power law density perturbations 
to an Einstein-de Sitter universe one expects a self-similar evolution. For the case of 
a Gaussian perturbation held the mean density prohle around high local maxima of the 
density perturbation held is given by the autocorrelation function. In the case of a scale 
free held whose power spectrum is given by P{k) oc where the relevant exponents in 
cosmology are in the range of —3 < n < 1, one hnds (HS): 

din) Qc ( 1 ) 

(Here is the initial radius dehned at some arbitrary time, well within the linear regime.) 
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A fiducial density profile is defined as the structure one would obtain by freezing all shells 
at their turn-around radius. This yield (HS): 

po(r)ocr“^°, (2) 


where 


7o 


3(3+ n) 
4 -I- n 


( 3 ) 


Assuming that the hnal (time averaged) radius r of the shell labeled by scales with r^, 
one hnds that the hnal density run, p(r) scales with the hducial one, namely p(r) oc 

In the case of a centralized perturbation, whose inner part is denser than its outer 


regions, one expects a secondary infall, namely the inner denser part collapses and relaxes 
hrst and the outer shells are collapsing on a longer time scale. As was pointed hrst by G77, 
the dynamics of the infalling less bound shells admits the product rM (r) as an adiabatic 
invariant, where M(r) is the mass that is momentarily interior to the radius r. Expressing 
the scaling of the hnal radius with the initial one by relating r/ to the turn around radius tq 
we write r/ = Ftq, where for a self similar infall F is a constant that depends on 79 . Now, 
the one-to-one dependence of tq on is easily found from energy considerations. Given 


all that, the adiabatic invariance implies that the rotation curve of the halos is given by: 


2 GM,(r,) GM„{r„) 

*rot „ I7'2„ 

Tf r 

Here M^iro) is the mass enclosed by the shell labeled at the turn around epoch. To 
complete this analytical calculation one should hnd the collapse, or contraction, factor F. 

For spherical systems, the mass interior to a given, i.e., Lagrangian shell, is conserved 
up to its turn around. However, beyond this stage shell crossing occurs and this mass is 
not conserved any longer. The equation of motion of spherical shells admits a self-similar 


solution which transforms the partial differential equation into an ordinary one, which 
is solved numerically (EG, Bertschinger 1985). These numerical exact solutions can be 
obtained analytically in the asymptotic limit of a shell which has past its turn around phase 
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much earlier than some given time, dehned as the ‘present’ time, and whose apocenter is 
much smaller than the current turn around radius, Rq. These asymptotic solutions have 
been solved subject to the approximations of spherical symmetry, scale invariance and the 
validity of adiabatic invariance (G77, FG, Bertschinger and Watts 1988). This problem 
has been recently addressed and an exact analytical solution for the collapse factor F 
in the asymptotic limit has been obtained (Zaroubi and Hoffman 1993, hereafter ZH). 
This solution is presented here and used in the comparison of analytical and numerical 
calculations. Gonsider a shell whose current apocenter (which is approximately equal to 
the time-averaged radius) is r and whose turn around radius is ro, then the collapse factor 
is given by: 


F 


r-C/ 


-1 


1 + (3 - 7o) 


'^°P(—)du 
u 


( 5 ) 


Here P{r jr') is the fraction of time a particle with apocenter r' spends inside r,U = Rq/to, 
Rq is the outer radius of the system or the present epoch turn-around radius and ro is 
the turn-around radius of the shell whose current apocenter is r. The crucial step is the 
evaluation of P{u). Various approximations were taken by G77, FG and Bertschinger 
(1985) and it has been solved exactly by ZH. Now, the integral in the denominator of 
Eq. 5 shows a very different behavior for 70 smaller or larger than 2 (which corresponds 
to n = —1, in the limit of self similar systems). In the limit of 77 3 > 1 the integral is 
dominated by its lower limit for 70 > 2 and by its upper limit for 70 < 2. Thus for 70 > 2 
the dynamics of the shell labeled by r is affected mainly by nearby outer shells such that 
r'>r, while in the other case of 70 < 2 it is dominated by shells that have turned around 
only recently, Tq^Rq and Tq is the turn around radius corresponding to r'. It follows that 
in the above limit one hnds self-similarity for 70 > 2. In the other case of 70 < 2 the 
hnal density run is a power law of 2 . Note, that in the case of 70 > 2 the effect of shell 
crossing is a local one, as only nearby shells of contributes to the mass interior to 
r. Therefore, in such a case one might hope that the model can be applied over a hnite 
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range of r, even if globally the strnctnre does not follow a power law and self-similarity. 
For 7 o < 2 the main contribntion to M(r) comes from the onter shells that have tnrned 
aronnd only recently. This introdnces a very sensitive dependence on bonndary conditions. 
The appropriate npper limit for the integral of Eq. 5 is t/ ~ 3 (ZH). For snch a valne the 
dependence of F on 70 is htted by: 

F = 0.186 + 0.15670 + 0.0137^ + 0.0177;^ - 0.00457(^ + 0.00327^ ( 6 ) 

The analytical model presented here relies on three basic assnmptions, namely spher¬ 
ical symmetry, scale invariance and the validity of the adiabatic invariance. Yet, for actnal 
halos emerging ont of Ganssian random helds it is clear that at least the hrst two as¬ 
snmptions do not apply. The strnctnre aronnd typical peaks is neither spherical nor is it 
characterized by a power law (BBKS) and it is certainly hnite. The qnestion thns arises 
is to what extent can the simple model describe the complex dynamical evolntion. Now, 
from previons nnmerical simnlations we know that at least in a statistical sense the model 
does correctly yield the hnal strnctnre of bonnd halos (QSZ, Frenk et al. 1988; for detailed 
comparison see Hoffman, 1988). Thns, applying the model to ensemble averaged density 
prohle reprodnces the ensemble average taken over all the objects in the simnlations, the 
dynamics of each one of which is properly calcnlated by the N-body codes. 

III. Numerical Simulations 

We model here the collapse of bound structures seeded by local density maxima of 
the primordial perturbation held in a hat Friedma nn universe. The perturbation held is 
assumed to be a Gaussian random held which is dehned by its power spectrum. For a scale 
free process the power spectrum is given by P{k) oc k'^. Given a particular realization of the 
perturbation held, constrained to be centered on a local density maximum, the subsequent 
dynamical evolution is calculated by three diherent methods. 

First, the hnal virial conhguration is calculated by using the SIM. This is done by 
taking the spherically averaged h-held (centered on the peak), and applying Eq. 4 to it. 



Thus, the analytical model that has been derived for scale-free initial conditions (ZH), is 
applied here to a general initial density prohle which is not self-similar. The collapse factor 
F is calculated here for each given spherical shell of radius r by using the local logarithmic 
derivative of the initial density prohle. The self-similar solutions are applied here locally, 
as if the shell crossing process depends on the local properties of the density held. In the 
case where the actual density prohle yields a hducial density which follows an power 
law with 7 o > 2, the actual 70 thus obtained is used to calculate the collapse factor F. 
However, for 70 < 2 we take it to be 70 = 2. The dependence of F on 70 is calculated 
by ZH. This relation has one free parameter, namely the ratio of the largest turn-around 
radius (at a given time) to the turn-around radius of the shell under consideration. Based 
on the N-body simulations here we estimate this ratio to be t/ = 3 in average. 

Second, the hnal virialized halo is calculated by evolving the initial system in full 
N-Body simulations. These N-body simulations are carried out by using the Treecode 
algorithm (Barnes and Hut, 1986, and Hernquist 1987, 1988, 1990). Third, we introduce 
a new numerical model to trace the evolution of high density peaks under the assumption 
of spherically symmetric force and realistic initial conditions, rather than the spherically 
symmetric force and initial conditions as assumed in the SIM. This code is used to trace 
the evolution of a single peak in an otherwise Einstein-de Sitter universe. Such a code 
stands in between the full N-body calculation and the simplihed analytical model. It turns 
out that the rotation curves of the virialized systems calculated using this method or with 
the Treecode, agrees well over most of the halos, excluding the very inner regions (see §1Z). 
Therefore this method is used here as a reference for testing the SIM predictions for most 
of the simulations. 

The location of the center of the peak is very crucial for this kind of approach and 
it should be carefully dehned; here we chose it to coincide with the point of maximal 
density. The time step used in this code is taken to be a hxed fraction from the dynamical 
time (rdyn) of the system. Where dt = lQ~'^Tdyn gives a very good energy conservation 
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< 1%) over a complete simulation. A softening parameter is used to avoid an infinite 
force at the origin. The simulated region has been chosen to be a sphere of radius R = 1, 
which along with fixing the gravitational constant to unity (G = 1) define the physical 
units of the simulation. The softening parameter of most of the simulations is 0.1. 

The two basic limitations of the present simulations are the nature of boundary con¬ 
ditions and the rather small number of particles. Here vacuum boundary conditions are 
assumed and therefore tidal interactions are ignored. As for the number of particles, in 
most runs N = 2000 particles are used. For one of the runs this has been increased to 
N = 4000 and 8000, in such a way that the increase in the number of particles was used 
to increase the extent of the simulated objects and not the code resolution. Now, these 
limitations certainly affect the outcome of the simulations. However, given that the sim¬ 
ple model has already been successfully confirmed by large scale simulations [e.g., QSZ), 
in which boundary conditions are properly handled, and given that our simulations are 
aimed at studying in a controlled way the dynamics of individual objects, we think that 
the present calculations are adequate for addressing the problems stated here. 

IV. Initial Conditions 

The formalism of constrained realizations of Gaussian random fields (Hoffman & Ribak 
1991) has been used to set the initial conditions. Here we are interested in making a 
particular realization of a Gaussian perturbation field, which is constrained to have a 
3a density peak located at the center of the computational box. A peak is specified 
by 10 constraints, namely the peak amplitude, the three first partial derivatives which 
are set to zero, the three eigenvalues of the second-order derivatives which form a 3 x 3 
symmetric matrix, and the three angles which define the direction of the eigenvectors. The 
three eigenvalues which define the shape of the peak are taken to have their mean values 
(BBKS), and the eigenvectors are taken along the Gartesian axis of the computational box. 
The realizations are designed to have such peaks at the center, but they have also other 
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unconstrained local maxima and minima. 


The physical systems simulated here are self-similar in the mean. The background 
universe is Einstein-de Sitter and the power spectra considered here are scale free, P{k) oc 
{n = —2, —1, 0,1}. Two sets of random numbers, which specify the phases and ampli¬ 
tudes and hence the given realization, are constructed here and are used for the different 
power spectra considered. The two realization are labeled as a and b. The number of 
particles used for each simulation is typically ~ 2000 except in two cases where we used 
4000 and 8000 particles with the hybrid code to test the effect of larger systems. 

A summary of the numerical experiments calculated here is given in Table I. The 
simulations vary with respect to power spectrum (n), family a or 6 of realizations, num¬ 
ber of particles, and the kind of code used. All models were calculated with the hybrid 
(symmetrical) code, some of which were calculated also by the full Tree-N-body code. 

[ Table 1 ] 

V. RESULTS 

Here we are interested in the gross density structure of the forming halos and this is 
described by the rotation curves, thus ignoring much of the hne structure. The resulting 
rotation curves calculated in the three different ways are given in Figs, l(a-c) of three 
models of family b with power spectra of n = —2,—1,0. Other models for which the 
rotation curves are calculated by the (analytic) secondary infall and the hybrid numerical 
algorithms are given in Fig. 2(a-c). A close inspection of the calculated rotation curves 
shows, hrst, the good agreement between the full N-body and hybrid codes over most of 
the mass of the halos. There is a discrepancy at the very inner region of the systems. 
The comparison of the analytically expected and numerically calculated curves is more 
difficult to interpret. In cases where the local effective 70 is smaller than 2, the rotation 
curve is predicted to increase with radius, corresponding to a diverging mass prohle, and 
the structure is determined by the latest infalling particles. Not surprisingly the numerical 
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simulations cannot reproduce such rising rotation curves; this has also been conhrmed 
by higher resolution N-body experements of QSZ and Crone et a/.(1994). However, for 
models with spectral index larger than —1 the effective 70 increases and a better agreement 
is expected. Indeed, in the case of the two realizations (a, 6 ) of n = 0 and in particular 
n = 1 the predicted rotation curves are close to the simulated one. Yet, in all cases the 
numerically calculated curves fall below the analytical ones, indicating a larger collapse 
factor (F) than the calculated factor. 

The N-body simulations of halos with predicted rising rotation curves poses severe nu¬ 
merical problems because of its extreme sensitive dependence on the boundary conditions. 
Previous simulations of a representative computational box that includes many halos and 
an n < —2 power spectrum yielded non increasing rotation curves [e.g., QSZ). To study 
this problem the n = —2 model ( 6 ) has been further simulated with an increasing number 
of particles, N = 2000,4000 and 8000, using the hybrid code. The added particles are 
used to simulate a larger initial volume and not to increase the dynamical resolution. The 
outcome is presented in Fig. 3, where the calculated rotation curves are presented and 
compared with the expected curve for the N = 8000 simulation. The change of the simu¬ 
lated curves with N is expected. In the N = 8000 case we hnd a flat rotation curve over 
about the inner half mass that is followed by a decreasing curve. Note that the predicted 
curve is now dominated by the outer shells and is decreasing. This manifests the sensitive 
dependence on the boundary conditions and shows that the SIM fails in the 70 < 2 case. 
The model is validated, however, in the 70 > 2 regime of declining curves. 

A possible explanation of this surprising agreement between the analytic model and 
the simulations lies in the ordered behavior of the collapsing systems in energy space as 
opposed to the random behavior in real (conhguration) space (Hoffman, 1988). As was 
noticed by Quinn and Zurek (1988) the seemingly very complicated collapse process looks 
very ordered and ’gentle’ when viewed in energy space. This leads us to the key question 
to be addressed here, namely to what extent a collapsing system ’remembers’ its initial 
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conditions. In particnlar we focns here on the energy (per nnit mass) behavior of individnal 
particles, which changes in time as {d/di){(^-\-/2) = {d/dt)4>. Therefore in collapsing 
systems the energy of particles changes in time, however it is not clear a priori whether this 
change is coherent or leads to a snbstantial phase mixing. This is stndied here by drawing 
scatter plots of the hnal vs. the initial energies of individnal particles in all the hybrid-code 
simnlations (Figs. 4 (a-f)). The fnll Tree-N-body simnlations yield similar scatter plots. A 
strong correlation exists in all cases, indicating a coherent and ordered evolntion in energy 
space. It is clear that snch systems do not go throngh an efficient phase mixing, and thns 
they retain the memory of their initial conditions. Note that there is a clear co nn ection 
between this correlation and the spectral index n, the tightest correlation corresponds to 
the n = 1 case and it decreases with n. This trend coincides with the dependence of the 
agreement between the SIM and the n um erical simnlations on the spectral index. 

To fnrther stndy the energy evolntion we rank the particles according to their binding 
energies, from the least to the most bonnd particle. The possible change of this rank¬ 
ing thronghont the collapse and virialization process is stndied by plotting histograms of 
the freqnency of the relative change of this ranking of all particles of a given hybrid-code 
simnlation Figs. 5(a-g). One hnds that in all cases some 70% of particles do not change 
their rank by more than 10%, and a very similar resnlt is fonnd in the fnll Tree-N-body 
simnlations. The approximate conservation of the ranking thronghont the dynamical evo¬ 
lntion proves that the cosmological collapse proceeds rather ’gently’ and involves no violent 
relaxation. 

Next, the trajectories of individnal particles are considered. In Figs. 6(a-b) the 
trajectories in the (r — t) plane of 20 randomly selected particles from the (-l;b) model 
are presented. Note, that this is the most snccessfnl model in terms of reprodncing the 
nnmerically calcnlated rotation cnrve. Yet, these trajectories are in complete disagreement 
with the ones calcnlated by FG and envisaged by the simple analytic model. Very few 
trajectories do show the gradnal shrinking of the tnrn-aronnd radins. In the light of 
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this disagreement, the success of SIM in predicting the the correct rotation curves seems 
to be very surprising. However one can argue that the rotation curves depend mainly 
on the energy distribution, which is much more smoother and more rounded than the 
underlying density distribution, that can account for the agreement between the model 
and the simulations (see §VI for more detailed discussion). 

VI. DISCUSSION 

The SIM has been rigorously formulated for self similar systems, yet in the canonical 
cosmological model structure arises from a random perturbation held and proto-structure 
are necessarily hnite and not self-similar. The question that is then naturally asked is 
to what extent this highly symmetric model is applicable to ‘realistic’ initial conditions. 
Previous numerical work [e.g., QSZ) has already proved the validity of the model when 
applied to N-body simulations and that in general the simple model correctly predicts the 
structure of bound objects, at least in the statistical sense. (This holds in the n > —1 case.) 
Yet, closer inspection of these simulations shows that the actual collapse process seems to 
be in complete disagreement with the one predicted by the simple model. It is this ambiva¬ 
lent behavior that has been studied here. Using the algorithm of constrained realizations 
special care has been given to the setting of the initial conditions of the simulated objects. 
Given the initial conhgurations, these were evolved dynamically to virial equilibrium in 
three different ways, namely by ‘exact’ N-body simulations, the analytical model and the 
hybrid N-body ‘monopole term’ code. Thus, a given initial conhguration has been evolved 
by three different methods and the hnal results are compared. Our basic conclusion is that 
within the limitation of the numerical simulations the gross agreement between the model 
and the large scale N-body simulation is co nfir med at the level of individual objects, in 
cases where such agreement is expected. 

Our main motivation here is to perform controlled numerical ‘experiments’ of the 
collapse onto local density maxima, with an emphasis on the controlled aspect. This is 
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achieved here in two ways. One is the nse of constrained realizations for setting the initial 
conditions. The other is the nse of the hybrid monopole term code that stands between 
the highly symmetric and analytic SIM and the fnll N-body simnlation. The comparison of 
the ontcome of this code with that of the standard Tree N-body code allows ns to separate 
the dynamical effects from that of (the lack of) spatial symmetry. We hnd here that within 
the technical limitation of onr simnlations, namely the lack of tidal interactions, the two 
modes of simnlations are in a good agreement. This proves that mergers mechanism does 
not stand at the basis of the snccess of the SIM. 

The basic physical nnderstanding that emerges ont of the present collapse ‘experi¬ 
ments’ is that of two different dynamical regimes depending on the logarithmic derivative 
of the (spherically averaged) fractional overdensity as a fnnction of radins, i.e., the effec¬ 
tive power law index n or eqnivalently the effective 70 of the hdncial tnrnaronnd density 
prohle. In the cosmological context these regimes correspond to a primordial pertnrbation 
held whose power spectrnm is dominated by the high wave nnmber modes n > — 1 or low 
wave nnmber modes n < — 1. In the regime of highly correlated pertnrbation held, n < — 1 
{i.e., in the mean 70 < 2 ), there is a strong dependence on the bonndary conditions, the 
dynamics strongly depends on the last collapsing shells. In snch cases the energy of par¬ 
ticles changes violently in time and one expects to reach at the end a state of statistical 
eqnilibrinm, mnch in the spirit of the Violent Relaxation proposed by Lynden-Bell (1967). 
It is important to note that the two very diherent approaches, namely the violent relax¬ 
ation and the SIM, predict the same hnal strnctnre of an asymptotic density prohle in 
the limit of inhnite spherical system. In the other regime of n > —1 the primordial strnc- 
tnres are sharply peaked and the dynamics of already collapsed shells is hardly ahected by 
the ongoing collapse of more distant and less bonnd shells. In snch a case order is being 
preserved in energy space, the initial conditions are well ‘remembered’ by the system and 
the hnal strnctnre rehects the initial conditions in a manner described here. 

The regnlar and ordered evolntion in energy space snggests a possible explanation as 
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to why a model based on spherical symmetry predicts correctly the hnal mass distribution 
of collapsing halos. There is a high correlation between the initial and hnal energies of in¬ 
dividual particles, and in particular the ranking in energy space is roughly preserved. Now, 
the virial structure depends on, and is basically determined by, the hnal energy distribu¬ 
tion, and consequently it depends on the initial energy distribution. In the linear theory 
of gravitational instability the total energy is determined by the gravitational potential, 
which is always much smoother and more rounded than the underlying mass distribution. 
Thus, the simple analytic model which seems at hrst glance to (incorrectly) describe the 
mass distribution, provides quite a good description of the initial energy distribution of 
particles, which in turn determines the hnal structure, in the case of 70 < 2 . 

The results presented here, and in particular the close agreement with the prediction of 
a spherical top-hat model, seem to be in conhict with earlier calculations on the role of shear 
in the gravitational collapse. These show that in the quasi-linear regime of gravitational 
instability in an expanding universe the shear (in the velocity held) accelerates the collapse 
and acts as a source of gravity (Hohman 1986 and 1989, Zaroubi and Hohman 1993b, 
Bertschinger and Jain 1994, Eisenstein and Loeb 1995). The shear, which is primarily 
induced by the tidal held, represents a deviation from spherical symmetry and thus its 
ehect seems to be inconsistent with the present results. Based on this Bertschinger (1994) 
postulated that local density maxima are not the sites where collapse occurs hrst. This, 
and the numerical simulations of Katz et al. (1993), shed doubt on the (linear) peaks - 
halos association. Yet, the more recent simulations of van de Weygaert and Babul (1994), 
which were designed to study the ehect of shear, have reaffirmed that high enough {^2a) 
peaks evolve to form halos. They found, however, that shear ahects the outer envelopes 
of halos. The simulations presented here were not designed to study these ehects and the 
limited dynamical range does not allow proper modeling of the external shear, however the 
internal shear which arise from the aspherical local matter distribution is well presented. 
This affects the collapse signihcantly in the quasi linear regime (Zaroubi and Hoffman, 
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1993b), yet the final virial structure is consistent with the top-hat model prediction. Note 
also that in N-body simulations where the computational box contains a typical patch 
of the universe with many proto-objects, both internal and external shears are properly 
sampled. Yet, such simulations basically confirm the predictions of the simple, top-hat 
spherical model (QSZ, Frenk et al. 1988, Crone, Evrard and Richstone, 1994). A tentative 
conclusion is that shear affects the dynamical evolution in the quasi-linear regime only, 
but the virial structure is hardly affected by it. However, this should be further tested in 
controlled N-body simulations. 
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TABLE 1 


The Parameters Of The Simulations 

Label 

Power Spectrum 

Family 

Particles’ ^ 

Simulations 

0 ;a 

n= 0 

a 

2000 

Symm.* 

l;a 

n= 1 

a 

2000 

Symm. 

-2;b 

n=-2 

b 

2000 

Symm. + Nbodyf 

-l;b 

n=-l 

b 

2000 

Symm. + Nbody 

0 ;b 

n= 0 

b 

2000 

Symm. 

l;b 

n= 1 

b 

2000 

Symm. + Nbody 

-l;b,4 

n=-l 

b 

4000 

Symm. 

-l;b,8 

n=-l 

b 

8000 

Symm. 


t Nbody, indicates full N-body simulation 
* Symm., indicates hybrid code simulation 
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Figure Captions 


Figure 1. 


Figure 2. 


Figure 3. 


Figure 4. 


Figure 5. 


Figure 6. 


The final rotation cnrves calcnlated in the three methods, the theoretical model 
(dashed lines), the hybrid monople-term code (long dashed lines), and the fnll N- 
body code (solid lines). The cnrves are drawn for the models labeled (see table I) (a) 
-2;b, (b) -l;b, and (c) l;b 

The final rotation cnrves calcnlated theoretically (dashed lines) and nnmerically by 
the hybrid code (long dashed lines). The cnrves are drawn for the following models 
(a) 0;a, (b) 0;b, (c) l;a. 

Rotation curves of the n = —2 model (labeled -2;b) calculated with the hybrid code 
for N=2000 particles (solid line), N=4000 (dotted line), and N=8000 (small dashed 
line). This is compared with the theoretically expected rotation curve for N=8000 
(long-dashed line). Note that the length scale here differs from that of Fig. la. 

A scatter plot of the final energies compared with the initial ones of the bounded 
particles as calculated with the hybrid code. The different graphs corresponds to the 
following models (a) 0;a, (b) l;a, (c) -2;b, (d) -l;b, (e) 0;b and (f) l;b. The correlation 
coefficient is denoted on each graph. 

Histograms of the change in the ranking (by binding energy) of the particles in the 
energy space, as calculated by the hybrid code. The histograms show the percentage 
of particles which changed their ‘energy’ ranking within a certain percentage. The 
models here are the same as in Fig. 4. 

The radial trajectories of randomly selected particles tin the Tree-N-body simulation 
for the model l;b. All radii are normalized by the maximal turn-around radius of each 
model. 
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